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We present a microscopic model for Na clusters embedded in raregas matrices. The valence electrons of 
the Na cluster are described by time-dependent density-functional theory at the level of the local-density 
approximation (LDA). Particular attention is paid to the semi-classical picture in terms of Vlasov-LDA. The 
Na ions and Ar atoms are handled as classical particles whereby the Ar atoms carry two degrees of freedom, 
position and dipole polarization. The interaction between Na ions and electrons is mediated through local 
pseudo-potentials. The coupling to the Ar atoms is described by (long-range) polarization potentials and 
(short-range) repulsive cores. The ingredients are taken from elsewhere developed standards. A final fine- 
tuning is performed using the NaAr molecule as benchmark. The model is then applied to embedded 
systems NasArjv. By close comparison with quantum-mechanical results, we explore the capability of the 
Vlasov-LDA to describe such embedded clusters. We show that one can obtain a reasonable description by 
appropriate adjustments in the fine-tuning phase of the model. 
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1 Introduction 

Structure and dynamics of metal clusters has been a much studied topic over the last two decades, see e.g. 
[1, 2, 3]. While optical response was in the focus of numerous earlier studies, recent developments aim at 
exploring cluster dynamics in all its aspects [4]. The preferred theoretical tool for describing clusters in all 
dynamical regimes is the time-dependent local-density approximation (TDLDA) which treats the cluster 
electrons at a fully quantum mechanical level, however in an independent-particle picture as motion in a 
common self-consistent mean-field. For highly excited systems, one can employ a semi-classical approx- 
imation in terms of phase-space dynamics described by the Vlasov equation, in the cluster context called 
a Vlasov-LDA. The Vlasov equation was originally used in a genuinely classical environment for plasma 
physics [5]. Nuclear physicists have adapted it to finite quantum systems for the description of violent 
heavy-ion collisions [6]. The motivation is twofold: first, the Vlasov description becomes more efficient 
for highly excited system, and second, the semi-classical approach makes it feasible to include correla- 
tions simply in terms of an Uhling-Uhlenbeck collision term (extended Vlasov to VUU). The successful 
applications in nuclear dynamics have motivated the transfer of the method to cluster physics, in the first 
round using a jellium model for the ionic background [7, 8, 9] but later also in connection with detailed 
ions [10, 11]. These applications of Vlasov-LDA and VUU dealt up to now only with free metal clusters. 
Many experiments, however, are done on supported or embedded clusters, see for example the study of 
the systematics of optical properties in small Ag clusters [12]. There is thus a need to develop reliable and 
efficient models for the cluster-substrate interface in connection with cluster dynamics in various regimes. 
Such developments are underway for TDLDA in embedded clusters whereby raregas matrices are consid- 
ered as a substrate with simple constituents [13, 14, 15]. They take up the modeling of the Ar atoms and 
their interaction with the Na cluster from earlier quantum-chemical considerations [16]. It is the aim of 
the present paper to investigate the performance of Vlasov-LDA and VUU for metal clusters embedded in 
raregas matrices. To that end, we recycle the interface model from the TDLDA studies and combine it with 
the Vlasov-LDA description for the cluster electrons. We intend to develop a parameterization of the model 
which provides acceptable ground state structures. These are then used as laboratory for dynamical studies. 
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electrons 


{<p n {r),n=l...N cl } /(r,p) 


cluster ions < — ► 


{R / ,/ = 1...7V ion } 


Ar core 


{R ,a = 1...JVA,} 


Ar valence cloud < — ► 





Table 1 Short summary of the degrees of freedom of the model. Note that the center of the valence cloud R' a serves 
to characterize the Ar dipole moment through d = eqA r s = e<ZAr(R' — R) where qAr is the net charge of the valence 
cloud. 

In the first exploration here, we use a simple excitation mechanism (instantaneous boost) to evaluate the 
basic features of dynamical response in various regimes, namely spectral distributions and ionization. The 
validity and limitations of Vlasov-LDA are scrutinized by comparison with results from quantum mechan- 
ical TDLDA. We also exploit the option of VUU for a first exploration of dynamical correlations and their 
interlay with the substrate. 

The paper is outlined as follows: In section 2, we present the model in detail thereby providing all in- 
formation for full TDLDA as well as for the semi-classical approaches. In section 3, we summarize briefly 
the emerging equations of motion and their numerical solution. In section 4, we discuss the calibration of 
the model using the NaAr molecule as benchmark. Finally in section 5, we present and discuss first results 
for an embedded cluster with Na 8 ArAr as test cases. 

2 The Na-Ar energy functional in detail 

2. 1 The degrees of freedom 

The key constituents are the valence electrons of the metal cluster. Starting point for their description is 
quantum mechanical mean-field theory where each electron is associated with one single-particle wave- 
function {(p n (r),n = l...JV e i}. The tp n are determined by a (time-dependent) Kohn-Sham equation. 
The manifold of occupied wavefunctions is comprised in the one-electron density matrix p(r, r') = 
S n Pn( r )'Pn ( r ')- We perform a semi-classical approximation in order to obtain a description in terms 
of the one-electron phase-space distribution /(r, p) [17]. This yields the Vlasov-LDA which has shown 
to provide a reliable description of dynamical processes in free metal clusters with detailed ionic structure 
[18, 19]. The other ingredients are the Na ions in the cluster and the raregas atoms. They are treated by 
classical molecular dynamics in terms of positions R and associated momenta P. The Ar atom carries 
as internal degree of freedom a dipole moment to account for the polarization interaction. Effective Ar 
potentials for static calculations contain that implicitly. True dynamics requires explicit dipoles as we use 
it here. It is convenient to formulate the Ar dipoles d a in terms of the displacement R' a — R a = d a /eqAi- 
All these degrees of freedom are summarized in table 1 . 

The detailed description of the Ar atoms is outlined in the next subsection 2.2. What the Na+ ions is 
concerned, we neglect the small dipole polarizability of the Na + core and treat it merely as a monopole. 
The ions are taken as point charges in the interaction with the Ar atoms. Soft and local pseudo-potentials 
are used for the interaction with the cluster electrons [20], see eq. (5). The density of the electrons is given 
naturally as defined in mean-field theories. It reads in the quantum mechanical version 

Pcl (r) = ]T|^(r)| 2 (la) 

a 

and in Vlasov-LDA 

p el (r)= f rf 3 p/(r,p) . (lb) 
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2.2 Construction of soft dipole potentials 



A finite dipole momentum in Ar atom is generated from Coulomb fields of the other constituents or from 
the excitation mechanism. In turn, it produces a dipole field to the outside world. However, a pure dipole 
field has a huge singularity at the origin. This is unphysical because we are not having point dipoles in 
practice. With a closer look at the Ar atom we see that the dipole is generated from a deformation (more 
precisely, displacement) of the outer electron shell (3s-3p-shell) and the finite size of the source delivers 
regular Coulomb fields everywhere. We thus specify the dipole momentum of the Ar atom in more detail. 
The Ar atom is thought of consisting out of a valence electron cloud with center of mass and of the 
corresponding positively charged core placed at the atom position R a . (In fact, one should use the center- 
of-mass of core plus valence cloud but the recoil effect on the Ar core is negligible, thus neglected). Both, 
valence cloud and core have the same absolute value of charge q\ r delivering together a neutral atom. A 
small displacement of the valence electrons against the core produces a net dipole moment. The profile 
of valence cloud and core is described by the same Gaussian with width <TAr, yielding together the charge 
distribution of one Ar atom as 



PAr,a(r) 



e<?Ar 



.3/2. 



Ar 



cxp 



(r-RJ 



'Ar 



exp 



(r- 



-K) 2 



'Ar 



(2a) 



where the index a stands for Ar atom number a implying atom position R a and valence cloud R' a . The 
Coulomb potential of this charge density produces the polarization-part of the potential exerted from the 
Ar atoms, 



<°>)=^Ar 



erf(|r-RJ/<7 Ar ) erf (|r-R^|/a Ar ) 



Ir-R„ 



Ir-R' 



(2b) 



where g Ar is the effective charge of the Ar cloud and core (see below). The error function therein is defined 

as 



erf(r) 



2 r 

Jo 



dxe 



(2c) 



2.3 The total energy 



The model is fully specified by giving the total energy of the system. It is composed as 



-^total = -^Na,el + ^Najon + -Eion,cl +£?Ar + -^coupl + £"vdW 



(3a) 
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-E'Na.ion 
-^ion,el 



= E kin + ^Hartrcc.Pel) + ^xc(Pel) + / drV ext (T, t)p e \(r) 



E 

i 

E 

i 

E 



2M N . 



E 

KJ 



R/ R,; 



J dry Ps p(|r-R/|)p e i(r) 



P' 



2M Ar 



^ 9m 



2m Ar 



-fc Ar (Hi - Ra) 



coupl 



+ E [/ ^p Ar , (r)^>) + v|X C) (R a - RaO 

a<a' 

+ E / dr Kxt(>,f)PAr : a(r) , 
El^^^ + ^NaAr^I-Ra)" 

+ [ dr Pcl (r) Y, ^aT>) + W elAr (|r - R a |)" , 
e ^EM]^ (/ dvf a (v)p el (v)^ - Jdvf a (vfp, 



i(r) 



= V 



erf(|r-RJ/a Ar ) 



|r-RJ 

The electronic kinetic energy depends on the level of approach, i.e. 

£ki„ = ^|V^„| 2 

n 

in the quantum case and 



d 3 p^-f(r,p) 
2m c \ 



(3b) 
(3c) 

(3d) 



(3e) 



(3f) 
(3g) 
(3h) 

(3i) 



(3j) 



in the Vlasov case. The Van-der-Waals energy i?vdw is a correlation from the dipole excitation in the 
Ar atom coupled with a dipole excitation in the cluster. It requires, in principle, an energy denominator 
involving the basic excitation energy AEat in the Ar atom and the typical excitation energy in the Na 
cluster which is w lumic- We exploit the fact that WMie *C A£?Ar and eliminate the energy denominator in 
the spirit of a closure approximation. This amounts to compute the variance of the dipole operator in the 
cluster. It employs, in fact, an effective dipole operator f a , see eq. (3h), which corresponds to the dipole 
field from the smoothened Ar charge distributions. The dipole variance as expectation value within the 
many-body wavefunctions is, furthermore, simplified in terms of the local variance. This yields finally the 
approximation (3g) to the Van-der-Waals energy. Its computation is still very cumbersome. It becomes 
particularly involved in the semi-classical case. On the other hand, Vlasov-LDA is intended for the regime 
of high excitations where subtle details, as e.g. the Van-der-Waals force do not matter. Thus we skip i?vdw 
in Vlasov-LDA. 

It is obvious that the polarization potentials are handled in terms of effective smoothened pseudo- 
densities (2a) and their resulting Coulomb potential (2b). The repulsive core potentials require explicit 
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modeling. We take it mainly from literature. The contributions are 

^iAr(r) = e \ + e t ( U l} (4a) 
V A?Ar\ R ) = e 2 1.367 -10- 3 





/?Na Q!Na ^Na CNa,6 CNa,8 


quantum mechanical LDA 
Vlasov-LDA 


1.7624 a^ 1 1.815 a 334.85 Ha 52.5 Ha ag 1383 Ha a§ 
1.670 a" 1 5.6 a 242.32 Ha 225 Ha ag 4410 Ha ag 



Table 2 Parameters for the ion-Ar pseudo-potential in Na according to the form (4c). 

For the Ar-Ar core interaction we employ a Lennard-Jones type potential with parameters such that 
binding properties of bulk Ar are reproduced. The Na-Ar core potential has been fitted according to [21]. 
Note that the Na-Ar potential from [21] does contain also a long range part oc ctAr which accounts for the 
dipole polarization-potential. We describe that long range part explicitely and have to omit it here. We thus 
choose the form (4c). The pseudo-potential VL c iAr for the electron-Ar core repulsion has been modeled 
according to the proposal of [16]. Its parameters will be re-tuned as discussed below. 

2.4 The ion-electron pseudo-potential 

The pseudo-potential for the coupling (3d) between Na ions and electrons is taken in a soft and local form 
as [20] 

Vp sP (r) = -- 
r 

with the error function given in eq. (2c). The parameters were adjusted to the properties of Na atom and 
bulk. This serves to provide simultaneously correct bond lengths for Na molecule and clusters as well as 
as good description of optical response, all in the framework of quantum-mechanical (TD)LDA [20]. A 
semi-classical description of Na clusters does also work very well with a pseudo-potential in the form (5). 
But it requires a new adjustment of the parameters suited for the Vlasov approach [22]. Both version for 
the parameters are summarized in table 3. 

2.5 The parameters for the polarization potential 

The polarizability of the Ar atoms is described by the model of a rigid cloud of valence electrons which is 
oscillating against the remaining raregas ion. The parameters of the model are: 





Cl 


C2 


cti 


CT2 


quantum mechanical LDA 


-2.29151 


3.29151 


0.681 a 


1.163 a 


Vlasov-LDA 


-0.5 


1.5 


0.7778 a 


1.5556 ao 



Table 3 Parameters for the ion-electron pseudo-potential in Na according to the local form (5). Different sets are to 
be used for quantum mechanical LDA and for semi-classical Vlasov-LDA. 

Copyright line will be provided by the publisher 



^6.501^ 12 ^6.501^ 6 



R 1 + e<W-R { R6 + R8 ) 



Cierf (7fc) +C2erf (vfc). 



6 



F. Fehrer 1 , M. Mundt 12 , P.-G. Reinhard 1,2 , and E. Suraud 2 : Modeling Na clusters in Ar matrices 



binding distance do 

depth E Q of the X 2 T, + ground-state 

energy of the transition X 2 £+^B 2 £+ 




Table 4 Basic properties of 
the NaAr molecule as they are 
used for fine tuning the model 
parameters. 

qAi = effective charge of valence cloud 

WAr = effective mass of valence cloud 

kAi = restoring force for dipoles 
These parameters are adjusted to reproduce the basic response properties of the raregas atom in the low 
energy domain, i.e. we choose to reproduce the static polarizability ana an d the second derivative of the 
dynamical polarizability d 2 ao |c^o- To that end, we make a one-pole model having the dispersion relation 



a Ar (0) 

-L 71 



. "Ar 2 
"Ar H 2"^ 



(6a) 



The frequency, in turn, is related to spring constant and mass as usual 



uj = ^k Al /m Al . 

On the other hand, the spring constant is directly related to the static polarizability as 



p 2 n 2 
e ffAr 

C*Ar 



We identify the mass of the electron cloud with its charge and the electron mass as 

flAr = qArm pA . 

thus we find for the effective charge of the valence cloud 



<7Ar = 



OArm-d^o 



(6b) 



(6c) 



(6d) 



The folding width of the valence cloud is determined such that the restoring force from the folded Coulomb 
(for small displacements) reproduces the spring constant &Ar- This yields 



org 

For Ar we use 



«A 



47T 



r 3(2^) 3 / 2 



1/3 



a Ar (0) = 11.08 a 3 , 



w 



1.755 Ry 



(6e) 



(6f) 



The eqs. (6) determine the parameters for the polarization potentials at the side of the Ar atoms. 



2.6 Fine tuning of electron- Ar core potential 

The parameters of the electron- Ar core potential (4a) determine sensitively the binding properties of Na to 
the Ar atoms. We use that potential as a means for a final fine-tuning of the model. The benchmark for 
adjustment is provided by the Na-Ar dimer. The data are taken from [23] and [24] and summarized in table 
4. There are three data points and three parameters to be adjusted, namely A c \, (3 C \, and r c \. It turns out 
that these three parameters compensate each other to some extent. One can, e.g. vary the potential height 
^4 e i and obtain about the same fit with re-tuning the other two parameters. Fortunately, we find a reasonable 
fit to the wanted data in spite of that restricted flexibility. Moreover, we exploit the freedom in the choice 
of A c \ to produce the softest reasonable core potential. It is obvious that such a subtle adjustment depends 
crucially on the level of approach. We thus obtain two different sets for quantum calculation and for Vlasov 
LDA, similar as it was the case already for the Na ion-electron pseudo-potential. The final parameters are 
summarized in table 5. 
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A e \ /3 e i r e i 


quantum mechanical LDA 
Vlasov-LDA 


0.47 1.6941 /a 2.2 a 
0.4 1.44 /a 1.8 a 



Table 5 Parameters for the ion- 
electron pseudo-potential in Na 
according to the local form (5). 
Different sets are to be used for 



quantum mechanical LDA and 
for semi-classical Vlasov-LDA. 



3 Dynamical equations and their solution 

3 . 1 The equations of motion 

The energy functional once specified (see section 2) determines structure and dynamics of the coupled 
system. This follows standard techniques throughout. We thus give here only a very brief summary. More 
details can be found in [4, 4] and specifically for the semi-classical approach in [9, 18]. 
The equations of motion are determined variationally. They read 

ldt<Pn = + U KS*j <fn , (7a) 

U KS (r) = , (7b) 

d t R = V P £ to tai , dtP = -V R £ t otai , R e {R /; R a ,Rj} , (7c) 

for the case of quantum-mechanical propagation of electrons. Only eq. (7a) is to be replaced by the Vlasov 
equation for the phase-space distribution /(r, p, t), i.e. 

ft/ + {^ + tw}=0 (7d) 

where {..,..} is the classical Poisson bracket. Note that the self-consistent Kohn-Sham potential U K s using 
the LDA energy functional is employed. That is why the scheme is called Vlasov-LDA. In fact, the above 
equations of motion describe electronic dynamics, in terms of TDLDA or Vlasov-LDA, coupled to classical 
MD for ions and atoms. The whole scheme is then TDLDA-MD in the quantum case or Vlasov-LDA-MD 
in the semi-classical approximation. 

It is the great advantage of the semi-classical approximation that it allows to include dynamical elec- 
tronic correlations through an Uhling-Uhlenbeck collision term. This extends Vlasov-LDA to the Vlasov- 
Uhling-Uhlenbeck scheme (VUU). The eq. (7d) is then extended to 

a t /+|^ + t/ KS ,/} =/uu , (7e) 



Juu (r, pi, t) = J d 3 P2 d 3 P3 d 3 P4 W(12, 34) - f™f°f) , (7f) 

W(12,M) = -\^6(p 1 +p 2 -p 3 -p 4 )5(E 1 +E 2 -E 3 -E 4 ) , (7g) 
m~[ du) 

= /(r, P< ,t)/(r, Pj ,t) , (7h) 

fST = (l-(2nhf^^)(l-(2,nfl^^l) , (7i) 



where i = (r, p) is used as an abbreviation and da/dui is the in-medium electron-electron cross section. 
We use here the cross section which had been evaluated for Na clusters [9, 18]. 

The numerical solution of the coupled equations of motion employs different methods for the differ- 
ent ingredients [25]. The quantum-mechanical Kohn-Sham equation (7a) is solved by the time-splitting 
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method [26]. The classical MD (7c) for ions and atoms is propagated by the (velocity) Verlet algorithm 
[27]. The handling of the Vlasov equation (7d) and associated collision term (7f-7i) is more involved. It 
will be detailed in the next subsection 3.2. 

The static solution has to be provided before any dynamics can take off. The ionic/atomic configuration 
is found by a mix of simulated annealing and direct gradient towards the energy minimum. Standard 
methods for the stationary Kohn-Sham equation apply in case of the quantum mechanical treatment. Most 
involved is again the semi-classical case. A stable ground state distribution / (r, p) for the Vlasov equation 
is determined by the Thomas-Fermi method [9, 18]. It is to be noted that the density and mean field entering 
Thomas-Fermi has to be properly augmented with the test-particle folding (see next section) in order to 
produce consistent input for the Vlasov solution with the test-particle method. 



3.2 Test particle method 

The phase-space distribution is a function in six-dimensional phase space. Its representation on a (six- 
dimensional) grid is even today beyond the capacity of computers. One thus uses a representation in terms 
of test-particles. This technique had been imported for cluster physics from the well developed nuclear 
physics case [6]. The smooth one-body phase-space distribution / is represented by a swarm of pseudo- 
particles each one having a certain width, i.e. 

N N ™ 

/(r,P,t) = -r 7 £L ^.gr(r-r l (t)) ffp (p-p l (t)) , (8a) 
iv PP i=i 

r 2 

g x (x) = (2trt 2 )- 3 / 2 exp ( - ^ ) . x e {r,p} , (8b) 

H , (8c) 
1.7836 ao . (8d) 



The folding widths a r and a p are free parameters of the method. Relation (8c) is connected to the 
Husimi picture [28] which gives the elementary distribution g r g p the interpretation of a minimal quantum- 
mechanical wavepacket being used as means to derive and justify the semi-classical limit in the smoothest 
way [17, 29, 30]. There remains ay as free parameter. It can be used to accommodate quantum-mechanical 
binding properties through the Vlasov approximation as good as possible and it has been used successfully 
that way for free clusters injellium approach [7, 8, 9] as well as with detailed ionic background [10, 11, 18]. 
We use the value (8d) which provides a good reproduction of free neutral Na clusters. 

Inserting the ansatz (8) into the Vlasov equation (7d) yields the following equation of motion for the 
test particles 

d t Ti(t) = Pi/m cl , (9a) 
dtPi® = -VnC/ KS , 9 , (9b) 
U K s, g = g r *U KS = J d 3 r'g r (r-r')U K s(r') ■ 



(9c) 



The Kohn-Sham potential is obtained from the ionic/atomic positions and the actual electron density com- 
puted with eq. (lb) which reads in the test-particle representation 

N Npp 

p el (r,t) = -^J2g r (r-r i (t)) . (9d) 

ly PP i=i 

Note that the g r folded Kohn-Sham potential (9c) is employed in the test-particle propagation (9b). This 
complies with the Husimi picture and it allows an efficient handling of density and mean-field on a spatial 
coordinate-space grid. The density is accumulated on the grid with eq. (9d) using the g r so to say as 
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interpolation functions. The Coulomb and exchange-correlation potentials are determined on the grid 
composing the Kohn-Sham potential together with the potentials from ions and atoms. The forces on the 
test-particle are retrieved from the grid again with the help of the folding functions g r . This is the practical 
side of eq. (9c). 

The finite representation in terms of test-particles induces a principle long-time instability towards a 
final Boltzmann equilibrium distribution. This is unphysical because it violates the Pauli principle [31,32]. 
The critical time depends on folding width and number of test particles. The folding width is fixed by the 
adjustment to quantum features. It remains the particle number. We chose N tp = 4000 per electron. 
This provides sufficient stability over the time scale considered (200 fs). The collision term (7f) in VUU 
contains the Pauli blocking. This provides long time stability of the propagation. 

The practical evaluation of the collision term is already simplified by the test particle. The nine-fold 
momentum-space integration is replaced by four-fold summation over test-particles, each factor /, con- 
tributing one sum. The scattering process connects two incoming momenta with two outgoing ones, 
(Pi)P2) — > (P3jP4)- The outgoing momenta are fixed by momentum and energy conservation laws 
up to a scattering angle. The scattering angle is chosen from the unit sphere at random. 

3.3 Excitation mechanism and observables 

Typical cluster excitations are promoted by laser pulses or collision with highly charged ions. Both can 
be described as external time-dependent Coulomb fields acting on the cluster. The coupling to the field 
Vext ( r j t) is provided for that purpose in the basic energy functional (3). A laser pulse, e.g., is characterized 
by V cv ± oc D cos (cut) /profile (t). In this first survey, we are interested mainly on general spectral properties. 
These are explored in the limiting case of an arbitrarily short laser pulse. It is realized by an instantaneous 
initial boost of the whole electron cloud which reads, in the Vlasov case, 

/(r,p,t = 0) = / (r,p + b) (10) 

where f is the ground state distribution and b is the boost momentum. In the quantum case the boost is 
applied to each electron orbital as 

if a (r) -^<p a {r,t = 0) =e'V(r) . (11) 

The size of b determines the dynamical regime. Small b explore the regime of linear response. 

The dipole spectrum (optical absorption) is obtained from the subsequent dipole signal D(i) = J d 3 r 
/9 e i(r, t) er. It is recorded over the whole time evolution. The dipole strength is finally obtained from the 
imaginary part of the Fourier transform, 9{D(cj)}. The three vector components correspond to the modes 
in x-, y-, and z-direction. For details see [4, 25, 33]. 

A further key observable is electron emission. The computation of the number of emitted electrons 
N esc is straightforward. In TDLDA, we use absorbing boundary conditions [25]. The lost electrons are 
determined as the complement of the still remaining electrons A osc = A(0) — J d 3 r p e i(r, t). In Vlasov- 
LDA, we draw a large sphere around the simulation area with radius 20 a and count explicitely all test- 
particles leaving that sphere. 

4 Benchmark NaAr 

Most of the parameters entering our model are taken from other grounds, i.e. from values published else- 
where and adjusted to independent data. The various constituents contribute strongly counteracting forces. 
This calls for some freedom to final fine-tuning. We exploit the anyway not so easily determined electron- 
Ar core potential for that purpose. The form (4a) is taken from quantum-chemical investigations [16]. 
The parameters are given free for fine-tuning. Thereby we consider the Na-Ar dimer as benchmark. The 
fully quantum mechanical TDLDA includes the Van-der-Waals interaction (3g) while the semi-classical 
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Fig. 1 The Born-Oppenheimer 
surface for the NaAr dimer 
computed semi-classically 
with the test-particle repre- 
sentation and compared with 
full TDLDA. The parameters 
for both TDLDA and semi 
classics are given in table 5. 
In the semi-classical case, two 
different core heights in the 
Ar-electron core potential (4a) 
are considered. The value 
A c = 0.4 is the final standard. 
The filled rhombus indicates the 
fitting point as listed in table 4. 



approach does not have that term. The missing contribution is then mimicked to some extent by the read- 
justment of the free parameters of the electron-Ar potential. Figure 1 shows the Born-Oppenheimer 
surfaces for the Na-Ar molecule. The TDLDA fit complies nicely with the data point. The semi-classical 
results manage to reach that approximately. The results from two different core heights serve to demon- 
strate the variances one can explore in the adjustment: while reducing the core height, one can shift the 
minimum closer to the wanted value, but at the same time one finds the bonding energy moving off the 
optimum. Whatever one does, one finds that the semi-classical Born-Oppenheimer surfaces fall off much 
more quickly than the TDLDA one. Note that the asymptotic in case of TDLDA is dominated by the Van- 
der-Waals attraction oc r~ 6 . Vlasov-LDA does not have that term and thus is bound to display a different 
asymptotics. Altogether, one sees that it is by no means trivial that a good adjustment can be found at 
all. Thus the TDLDA fit is a remarkable agreement and we are also satisfied with the semi-classical results 
which manage to find a fair agreement for the binding point. One has to keep in mind that the stronghold of 
Vlasov-LDA and VUU are the energetic processes. For these it suffices to provide a stable starting ground- 
state configuration which reproduces roughly realistic binding properties. This is achieved by the choice of 
parameters as listed in table 5. Thereby we have preferred the larger core height which provides the correct 
bond length while giving a slightly too small binding energy. Note that the semi-classical adjustment was 
done without explicit Van-der-Waals energy. This term is derived from quantum mechanical perturbation 
theory. Its implementation in a semi-classical environment requires deeper theoretical development and 
what we have found would increase the computational cost by an order of magnitude. That is not very 
useful for an approach which concentrates on processes far above the low energy scale of binding. 

A further benchmark point for TDLDA is the energy for the X 2 Y, + — ► B 2 T, + transition at the bonding 
configuration. The goal is a transition energy of 0.155 Ry (see table 5); the TDLDA result with the 
optimized parameters is 0. 164 Ry, a satisfying agreement. This benchmark cannot be applied for the semi- 
classical approach because there are no quantized single-particle states and thus no specific transitional 
energy spectra. We will apply Vlasov-LDA to metal clusters. These have the Mie surface plasmon as 
prominent excitation mode and this mode is determined by collective flow rather than by single-particle 
excitations. Thus we can expect a good reproduction of the Mie plasmon in spite of the missing single- 
particle adjustment. This is indeed the case as many previous applications for free metal clusters have 
exemplified [7, 8, 9, 10, 18, 11]. We have yet to see how the model performs for embedded clusters. 
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Fig. 2 Ionic and electronic 
radii of a Nag cluster embed- 
ded in Ar matrices of differ- 
ent size. Results from Vlasov 
and from quantum mechani- 
cal treatment are compared. 



5 The NagAr^v system 

5.1 Structure 

The ground state configurations involve 8-220 ions and atoms. This is hard to visualize in all detail. The 
basic structure, however, is simple [13, 14]. There exist three close isomers for the Nas clusters [14]. We 
consider here the Nag cluster in the shape consisting out of two rings with each four electrons. The Ar 
matrix was originally arranged in radial shells. It is slightly perturbed by the embedded cluster. But the 
shells still remain a useful guide and sorting scheme. Thus we will characterize the Nag cluster by its 
global root-mean-square radius and the Ar matrix is visualized as distribution of atoms in radial shells. 

Figure 2 shows the radii for the (embedded) Na 8 . The ionic radii are generally smaller than the electronic 
ones because the electronic cloud reaches farther out than the well localized ions. The difference between 
electronic and ionic radius is caused mainly by the Gaussian folding with width oy. It is thus nearly 
constant. It is of the same order (although somewhat smaller) as the difference for the quantum mechanical 
case. The latter shows an interesting detail: it is practically constant for cases with Ar matrix but jumps 
by 0.2 a for the free cluster. This is due to the onset of core repulsion, similar as the first step for Vlasov 
LDA. In the further evolution, the radii show a different trend. While they slowly grow with A^r in the 
quantum mechanical case, they shrink for Vlasov-LDA. This hints at a different mix of interactions with 
the Ar atoms. The Ar core repulsion is felt more strongly for the test-particle in Vlasov-LDA because 
these are associated rigidly with a Gaussian distribution. The repulsive core wants to push away the tail 
and thus pushes the whole Gaussian. This works much different for the quantum mechanical case. The 
electron cloud can react more "elastically". It allows for a dip near the Ar core and evolves else-wise almost 
unperturbed around it. Thus the dipole attraction can act more freely to pull the electronic tail outward. 

Figure 3 shows the Ar part of the combined systems in terms of the radial atom distribution. The 
differences between quantum LDA and Vlasov-LDA structures are dramatic for the small Ar matrices. 
They converge to each other with increasing system size, similarly to the convergence observed in the 
case of Na properties (Figure 2). The difference is again due to the rigidity of the fixed Gaussian electron 
distribution for Vlasov-LDA. This exaggerates core repulsion and diminishes dipole attraction. Moreover, 
we are missing the long range Van-der-Waals attraction (see figure 1). The pressure from the outside 
shells and the increasing dipole attraction from the far shells produce then the compression towards the 
quantum mechanical (and thus more realistic) structures. The example demonstrates the limitations of the 
test-particle method. The smooth Gaussian folding implied in the test-particle representation was perfectly 
suited for the soft ionic pseudo-potentials of a free metal cluster (and the more so for the jellium approach 
to a free cluster). But the inherent rigidity of the electronic test-particle distribution becomes a hindrance in 
connection with steep potentials as they appear here in terms of the repulsive Ar cores. The Vlasov-LDA is 
certainly not appropriate for detailed studies of subtle binding effects. We find, however, that the structures 
for larger matrices come out sufficiently well to provide a laboratory for dynamical studies in the non-linear 
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Fig. 3 The distribution of radial shell in Ar matrices of different sizes for Vlasov-LDA and for a fully quantum 
mechanical treatment. 



domain 


excitation energy 


linear 

semi-linear 

non-linear 


0.3 eV 
4.1 eV 
13.5 eV 



Table 6 Excita- 
tion energies for 
different domains. 



regime. Remember that VUU versus Vlasov-LDA is presently the only way to estimate microscopically 
the effect of dynamical electron correlations on the resonance excitations. 



5.2 Dynamics 

In order to explore the basic dynamical properties, we employ the most simple excitation mechanism of 
an instantaneous boost, as explained in section 3.3. It has only one parameter, the boost strength, which 
allows to scan the various dynamical regimes. We pick a typical strength for each regime as listed in 
table 6 and characterized by the initial excitation energy. To put that in perspective, one may compare it 
with the typical energies in the system: at the side of the Na cluster we have the ionization potential with 
IP« 3.5 - 4.5 eV and the energy of the Mie surface plasmon with w ( Mi °) w 2.6 - 2.9 eV. At the side 
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Fig. 4 Results of Vlasov cal- 
culations for the power spec- 
trum of the dipole mode along z 
axis for Na 8 (upper panel) and 
for Na8@Ari64 (lower panel) 
at different excitation energies. 
The spectra are normalized to 
peak maximum at 1 for better 
comparison. The dipole mo- 
ment has been recorded for 100 
fs. 

of the Ar atoms, we have their first electronic excitation and the IP above 20 eV. The linear regime is far 
below anyone of these scales, the semi-linear just of the same order, and the non-linear regime safely above 
(multi-plasmon, strong emission). 

The spectral properties from Vlasov-LDA (and VUU) are shown in figure 4 in terms of the power spec- 
trum of the dipole mode along z-axis (which corresponds to the approximate symmetry axis of the Na 8 
cluster). The general pattern are much similar to previous studies on Na<^ [34, 18]. There is a clear reso- 
nance peak under any conditions. It is not destroyed by high excitations. These well developed resonance 
features persist very well also for the embedded cluster. One may even find the pattern somewhat cleaner 
than for the free cluster. At second glance, one sees trends and developments. There is, e.g., the much 
larger width in case of VUU which is obviously generated by the electron-electron collisions. We will now 
work out the trends in terms of the key features of these simple spectra: the peak position cj( Mic ) and the 
width as full width at half maximum (FWHM). 

Figure 5 shows these basic features in case of the system Na 8 Ar„ in < 212) for quantum mechanical 
calculations, Vlasov-LDA, and VUU. The peak frequencies (two lower panels) depend only very weakly 
on the Ar matrix in all cases. The remaining small trend shows a red shift with increasing matrix size in 
the quantum mechanical results. That is due to the fact that the additional Ar shells enhance the attractive 
dipole interaction. The Vlasov results indicate the opposite trend. That is due to the decreasing radius of 
the inner Ar shell (see figure 3) which gives more weight to the increasing core repulsion. The changes with 
increasing excitation energy differ. The quantum mechanical results show a sizeable and systematically 
increasing blue-shift. It can be explained by the increasing ionization which, in turn, deepens the binding 
potential and thus increases the restoring force for the dipole oscillations [33]. The same trend is hinted 
in the Vlasov results for the free Na 8 shown in the upper panel of figure 4. But the Vlasov results for 
the embedded cluster (middle panel of figure 5) behave differently in that the step from the linear to the 
semi-linear regime shows a tiny (or no) red shift instead of a blue shift. The reasons for that are not yet 
fully clear. The step into the non-linear regime then produces the blue-shift. Finally, the step to VUU has 
negligible effect on the peak position as one would have expected. 

The uppermost panel of figure 5 shows the widths of the resonance peaks in the Vlasov cases. The 
spectral analysis uses a recording time of 100 fs and employs a filtering to avoid artefacts from insufficiently 
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plasmon (lower two panels) 
and its full width at half 
maximum (FWHM, uppermost 
panel) for Nas Ar n at various ex- 
citations as indicated. Lowest 
panel: quantum mechanical re- 
sults. Upper two panels: results 
from Vlasov-LDA (or VUU, re- 
spectively). 



relaxed signals [35]. This sets the energy resolution of the analysis to 0.1 eV. The width of the Vlasov 
results is hidden in the resolution except for the non-linear case which shows already some enhancement 
due to a dynamical spreading (through ongoing ionization) of the spectra. However, the collision term 
in VUU adds a substantial bit of 0.2 eV for free Na 8 and 0.1 eV for the embedded clusters. It is a bit 
surprising to see a reduction of the width when the matrix is around. This hints at a somewhat more 
harmonic potential for the embedded cluster. The effect deserves a deeper analysis which, however, goes 
beyond scope and limits of this paper. 

Figure 6 shows the time evolution of ionization (number of emitted electrons) in case of the non-linear 
excitation. The pattern for quantum-mechanical TDLDA and Vlasov-LDA are similar: there is a steep ini- 
tial increase from direct electron emission and then a quick bending over to an almost flat trend. The results 
agree even quantitatively for the free cluster. For the embedded system, on the other hand, the agreement 
may be just acceptable but is certainly less perfect with Vlasov-LDA showing about half of emission from 
TDLDA. The reason is that cores of the surrounding Ar atoms build up a narrow potential barrier of about 
1 eV height. This barrier reflects in Vlasov-LDA all test particles which do not have sufficient kinetic 
energy while quantum mechanical electrons in TDLDA can easily tunnel through. More violent excita- 
tions produce a larger fraction of faster electrons and will thus diminish that difference. This confirms 
once again the rule that semi-classical approaches become increasingly valid with increasing excitation 
energy. The VUU results in figure 6 show a qualitatively different trend. There remains a slope which 
hints at significant delayed emission of thermalized electrons. The free Na 8 shows the well understood 
deal: the VUU curve remains below Vlasov-LDA at early times because direct emission is suppressed in 
favor of internal thermalization while it crosses the Vlasov line at later times when thermal emission takes 
over [9]. In case of the embedded cluster, one finds the VUU emission enhanced in all time domains. We 
see here again a classical effect. The part of thermal emission in VUU is isotropically distributed. The 
electrons have thus more chances to escape through the saddles between the cores. Moreover, the more 
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Fig. 6 Time evolution of the num- 
ber of emitted electrons N esc af- 
ter strong initial boost with exci- 
tation energy 13.5 eV (non-linear 
regime). Results from quantum me- 
chanical TDLDA, Vlasov-LDA, and 
VUU are compared. Upper panel: 
for the system Nag@Ari64. Lower 
panel: for free Nas. 
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Fig. 7 Final number of emitted 
electrons N CBC after initial boost in 
semi-linear and non-linear regime 
(see table 6). Results from quantum 
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effective thermalization in VUU allows to fill more systematically the Boltzmann tail of the kinetic energy 
distribution until some electrons have gathered sufficient energy to surmount the extra barrier. 

Figure 7 shows a summary of asymptotic values for ionization for different regimes, methods and system 
sizes. There is little dependence on the size of the Ar matrix, except for the too small A^r = 30. Thus we 
recognize many features from the previous figure. We see again the strong suppression of (direct) emission 
in Vlasov-LDA. It is even more dramatic in semi linear regime (factor 8 suppression!) because it has a 
larger fraction of slow electrons which cannot surmount the first Ar barrier. We also find again a systematic 
trend of enhanced emission in VUU due to its larger contribution from isotropic emission and from the 
Boltzmann tail. Altogether, the results for electron emission hint that the propagation of slow electrons 
in the Ar matrix is much different between quantum TDLDA and semi-classical approaches. The latter 
still provide useful guidelines and presently the only chance to study effects of electron-electron collisions. 
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One should, however, take care to remain in the regime of strong excitations where slow electrons are less 
important. 

6 Conclusion 

In this paper, we have presented a model for the description of the dynamics of an embedded metal cluster 
(Na) in a rare gas matrix (Ar). The model employs a microscopic description of three basic ingredients: 
cluster valence electrons, cluster ionic cores, and raregas atoms. The cluster electrons are treated within 
time-dependent density functional theory (DFT). The cluster ions are treated classically. The raregas atoms 
are also treated as classical particles together with their dynamical polarizability as intrinsic degree of free- 
dom. This model allows to handle a large variety of dynamical situations at the side of the cluster, including 
electron emission when non linear excitations are considered. The cluster electrons can be treated at 2 lev- 
els of approximation within the realm of DFT: fully quantum mechanical time-dependent local-density 
approximation (TDLDA) and (semi-classical) Vlasov-LDA. It turns out that both approaches give results 
which are qualitatively, sometimes even quantitatively, comparable. This holds in the stationary as well 
as in the dynamical regime. Fully fledged TDLDA reproduces quite well experimental data and previous 
quantum chemistry approaches. The semi-classical approach also provides an acceptable description for 
larger raregas matrices, which is to some extend a surprise in view of the subtle energetic balance in these 
systems. 

The semi-classical approach becomes particularly suited when non-linear excitations are considered. It 
presently provides the simplest and practically only way to deal with dynamical correlations of the clus- 
ter electrons. This is achieved by complementing the original Vlasov equation by a Boltzmann-Uhling- 
Uhlenbeck collision term yielding the Vlasov-Uhling-Uhlenbeck (VUU) model. The collisional correla- 
tions induce significant differences between the Vlasov and VUU approaches when a sufficiently large 
energy is deposited in the system, for free clusters as well as for embedded ones. The test cases considered 
make it desirable to account for these dynamical correlations also in a quantum mechanical context. This 
is, however, an extremely difficult task still under development. In the mean time the VUU approach has 
proven to provide a valuable tool for investigations of the non linear dynamics of simple metal clusters 
embedded in a rare gas matrix. Further explorations of the various competing phenomena are underway. 
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